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Abstract 

This  report  deals  with  the  formulation,  implementation,  and  testing  of  three  numerical 
techniques  based  on  (i)  a  full  multiphase  approach,  (ii)  a  MUlti-SIze  Group  (MUSIG) 
approach,  and  (iii)  a  Heterogeneous  MUSIG  (H-MUSIG)  approach  for  the  prediction  of 
mixing  and  evaporation  of  liquid  droplets  injected  into  a  stream  of  air  flowing  inside  a 
combustion  chamber.  The  numerical  procedures  are  formulated  following  a  Eulerian 
approach,  within  a  pressure-based  fully  conservative  Finite  Volume  method  equally 
applicable  in  the  subsonic,  transonic,  and  supersonic  regimes,  for  the  discrete  and  continuous 
phases.  The  k-8  two-equation  turbulence  model  is  used  to  account  for  the  droplet  and  gas 
turbulence  with  modifications  to  account  for  compressibility  at  high  speeds.  For  the  purpose 
of  comparing  the  performance  of  the  various  methods,  three  configurations  involving  stream- 
wise  and  cross- stream  spraying  in  rectangular  and  cylindrical  domains  are  investigated  and 
solutions  for  evaporation  and  mixing  in  the  subsonic  and  supersonic  regimes  for  droplets 
sprayed  in  turbulent  flow  streams  are  generated.  Results,  displayed  in  the  form  of  droplet 
velocity  vectors,  contour  plots,  and  axial  profiles  indicate  that  solutions  obtained  by  the 
various  techniques  exhibit  similar  behavior.  Differences  in  values  are  relatively  small  with  the 
largest  being  associated  with  droplet  volume  fractions  and  vapor  mass  fraction  in  the  gas 
phase.  This  is  attributed  to  the  fact  that  with  MUSIG  and  H-MUSIG  no  droplet  diameter 
equation  is  solved  and  the  diameter  of  the  various  droplet  phases  are  held  constant,  as 
opposed  to  the  full  multiphase  approach. 
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Nomenclature 

Ap ) coefficients  in  the  discretized  equation  for  (J)'k  1 . 

Bp  1  source  term  in  the  discretized  equation  for  (p'k  1 . 

Br  breakup  rate. 

B'j  breakup  frequency. 

Cpk)  coefficient  equals  to  1 1  R(k)T(k) . 

CD  drag  coefficient. 

cP  specific  heat  at  constant  pressure. 

Cy  coalescence  rate. 

d  droplet  diameter. 

ds  Sauter  diameter. 

Dp  [0a>]  the  Matrix  D  operator. 

/  population  fraction. 

Fb  Body  force. 

FD  drag  force. 

h  static  enthalpy. 

h correction  correction  coefficient  for  heat  transport  in  droplet  model. 

Hp[(f){k  1  ]  the  H  operator. 

HP[u(t)]  the  vector  form  of  the  H  operator. 
correction  correction  coefficient  for  mass  transport  in  droplet  model. 


mass  rate  of  droplet  evaporation, 
volumetric  mass  rate  of  droplet  evaporation. 


n  number  density  distribution  function. 

Pk  production  term  in  k  and  e  equations. 

p  pressure. 

Pr  laminar  Prandtl  number  of  fluid/phase  k. 

Pi;  turbulent  Prandtl  number  of  fluid/phase  k. 

q  heat  flux. 

Q(k)  general  source  term  of  fluid/phase  k. 

R(k)  gas  constant  for  fluid/phase  k. 

Red  Reynolds  number  based  on  the  droplet  diameter. 

S  source  term. 

Sf  surface  vector. 

Sc  Schmidt  Number, 

t  time. 

T  temperature  of  fluid/phase  k. 

Uf  interface  flux  velocity  (v  'k  1  .S  / ) . 

u  velocity  vector. 

u,v  velocity  components  in  x-  and  y-direction,  respectively. 

Xjki  mass  fraction  due  to  coalescence  between  groups  j  and  k,  which  goes  into  group 
Y  vapor  mass  fraction. 

Greek  Symbols 

a  volume  fraction. 

[3(k)  thermal  expansion  coefficient  for  phase/fluid  k. 

p  Kolmogorov  micro-scale. 

5t  time  step, 

p  density. 

T  the  stress  tensor. 
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A  conductivity  coefficient. 

Ac  coalescence  efficiency. 

coc  collision  frequency, 

r  diffusion  coefficient. 

(p  general  scalar  quantity. 

A hv  latent  heat. 

A p  [(p(k)]  the  A  operator. 

Id,/liturb  laminar,  turbulent  and  effective  viscosity  of  fluid/phase  k. 

Q  cell  volume. 

oho2  turbulence  model  constants. 

Subscripts 

d  refers  to  the  droplet  discrete  liquid  phase. 

E  refers  to  energy  equation. 

eff  refers  to  effective  values. 

g  refers  to  the  gas  phase. 

i  refers  to  phase  i. 

k  refers  to  turbulent  kinetic  energy  equation. 

nb  refers  to  the  east,  west,  . . .  face  of  a  control  volume. 

NB  refers  to  the  East,  West,  . . .  neighbors  of  the  main  grid  point. 

P  refers  to  the  P  grid  point. 

s  refers  to  the  droplet  surface  condition. 
sat  refers  to  the  saturation  condition. 

£  refers  to  turbulent  eddy  dissipation  equation. 

co  refers  to  turbulent  eddy  frequency  equation. 

vap,g  refers  to  the  vapor  specie  in  the  gas  phase. 
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Introduction 

This  report  deals  with  the  formulation,  implementation,  and  testing  of  three  pressure-based 
numerical  techniques  for  the  prediction  of  mixing  and  evaporation  of  liquid  fuel  injected  into 
a  supersonic  stream  of  air  flowing  inside  a  combustion  chamber. 

Recently  there  has  been  a  revived  interest  in  the  injection  of  liquids  in  supersonic  streams, 
particularly  with  respect  to  fuel  injection  techniques  for  hypersonic  flights.  These  designs 
require  air-breathing  engines  capable  of  supersonic  combustion.  The  scramjet  (supersonic 
combustion  ramjet),  appears  at  present  to  be  a  practical  engine  for  these  types  of  applications 
since  it  is  capable  of  producing  useful  thrust  at  hypersonic  speeds,  using  supersonic  flow 
through  the  combustor.  The  scramjet  concept  itself  is  fairly  old,  and  was  the  subject  of  studies 
throughout  the  1960s  and  again  in  the  1980s.  However,  its  coming  to  fruition  depends  among 
other  things  on  the  development  of  numerical  tools  for  the  simulation  of  its  supersonic 
combustion  process  and  related  phenomena.  More  specifically  effective  penetration  and 
enhanced  mixing  of  hydrocarbon  fuels  in  a  gas  flowing  at  supersonic  speed  is  a  crucial 
ingredient  for  the  success  of  any  scramjet  design  [1],  Three  key  issues  govern  the 
performance  of  the  liquid  injection  process  in  the  scramjet  engine,  namely:  the  penetration  of 
the  fuel  into  the  free-stream,  the  atomization  of  the  injected  fuel  drops,  and  the  level  of 
fuel/air  mixing  [2],  It  is  important  for  the  fuel  to  penetrate  effectively  into  the  free-stream  so 
that  the  combustion  process  produces  an  even  temperature  distribution  otherwise  it  will 
mostly  occur  along  the  surface  of  the  combustor,  causing  inefficient  combustor  operation  and 
increased  cooling  problems.  Rapid  atomization  of  the  fuel  is  also  required  for  efficient 
combustion.  Increased  atomization  of  the  liquid  fuel  results  in  increased  fuel/air  mixing  which 
allows  a  higher  percentage  of  the  fuel  to  be  burnt  in  the  short  time  before  the  entire  mixture 
passes  out  of  the  combustor  (generally  the  flow  residence  time  is  of  the  order  of  few 
milliseconds  [3]).  This  paper  is  aimed  at  developing  three  numerical  methods  capable  of 
predicting  the  spreading  and  evaporation  of  liquid  droplets  injected  in  gases  flowing  at  all 
speeds. 
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The  complex  multi-phase  flow  phenomenon  governing  liquid  injection  applications  involves 
a  continuous  gas  phase  usually  composed  of  air  and  the  evaporating  vapor  species  from  the 
fuel  and  one  or  more  dispersed  liquid  phases.  Moreover,  the  phenomenon  also  entails 
turbulent  effects.  Approaches  for  the  simulation  of  droplet  transport  and  evaporation  in 
combustion  systems  can  be  classified  under  two  categories,  namely  the  Lagrangian  and 
Eulerian  methods.  Within  both  methods  the  gaseous  phase  is  calculated  by  solving  the 
Navier-Stokes  equations  with  a  standard  discretization  method  such  as  the  Finite  Volume 
Method. 

In  the  Lagrangian  approach  [4,5,6],  the  spray  is  represented  by  discrete  droplets  which  are 
advected  explicitly  through  the  computational  domain  while  accounting  for  evaporation  and 
other  phenomena.  Due  to  the  large  number  of  droplets  in  a  spray,  each  discrete  computational 
droplet  is  made  to  represent  a  number  of  physical  droplets  averaging  their  characteristics.  The 
equations  of  motion  of  each  droplet  are  a  set  of  ordinary  differential  equations  (ODE)  which 
are  solved  using  an  ODE  solver,  a  numerical  procedure  different  from  that  of  the  continuous 
phase.  To  account  for  the  interaction  between  the  gaseous  phase  and  the  spray,  several 
iterations  of  alternating  solutions  of  the  gaseous  phase  and  the  spray  have  to  be  conducted. 
Therefore,  the  computational  effort  for  strongly  interacting  two-phase  flows  with  the 
Lagrangian  method  is  rather  large.  Furthermore,  for  turbulent  flow  simulation  the  above 
model  has  to  be  augmented  with  a  stochastic  or  Monte-Carlo  approach. 

In  the  Eulerian  approach  [5, 6, 7, 8],  the  evaporating  spray  is  treated  as  an  interacting  and 
interpenetrating  continuum,  in  analogy  to  the  continuum  approach  of  single  phase  flows,  each 
phase  is  described  by  a  set  of  transport  equations  for  mass,  momentum  and  energy  extended 
by  interfacial  exchange  terms.  This  description  allows  the  gaseous  phase  and  the  spray  to  be 
discretized  by  the  same  method,  and  therefore  to  be  solved  by  the  same  numerical  procedure. 
Because  of  the  presence  of  multiple  phases  a  multiphase  algorithm  is  used  rather  than  a 
single-phase  one. 
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In  this  work,  the  numerical  foundations  for  the  simulation  of  supersonic  droplet  evaporation 
and  scattering  are  developed.  This  is  achieved  by  following  a  Eulerian-Eulerian  approach 
using: 

1-  the  full  multiphase  flow  model, 

2-  the  Multi-Size-Group  model  (MUSIG),  and 

3-  the  Heterogeneous  Multi-Size-Group  model  (H-MUSIG). 

The  three  approaches  are  implemented  within  an  all  speed  pressure-based  finite  volume  flow 
solver  in  which  a  droplet  evaporation  model  is  implemented.  The  k-e  [9]  two-equation 
turbulence  model  has  been  used  to  account  for  the  droplet-gas  turbulence  with  modifications 
for  supersonic  flows.  Droplet  turbulence  is  estimated  using  an  algebraic  model  based  on  the 
Boussinesq  approach  [7,10,11].  Different  droplet  sizes  are  considered  and  droplet  breakup 
and  coalescence  [12,13]  are  modeled  and  their  effects  incorporated  into  the  conservation 
equations  via  source  terms.  The  use  of  an  Eulerian  approach  has  many  advantages:  same 
validated  numerical  procedure  used  for  all  phases,  ease  of  implementation  of  acceleration 
techniques,  and  improvements  to  code  can  be  carried  over  to  all  phases. 

In  the  remainder  of  this  document,  the  governing  conservation  equations  for  both  gas  and 
liquid  droplet  phases  and  their  discretization  procedure  are  first  presented.  This  is  followed  by 
a  detailed  description  of  the  three  solution  methodologies.  The  resulting  algorithms  are  used 
to  solve  three  physical  configurations:  (i)  streamwise  injection  into  a  subsonic  and  supersonic 
stream  flowing  in  a  rectangular  domain,  (ii)  cross-streamwise  injection  into  a  supersonic 
stream  flowing  in  a  rectangular  domain  are  given,  and  (iii)  streamwise  injection  into  a 
supersonic  stream  flowing  in  a  cylindrical  domain.  Results  presented  in  the  form  of  droplet 
and  gas  velocity,  pressure,  gas  temperature,  vapor  mass  fraction,  air  volume  fraction,  and  gas 
turbulent  viscosity  fields  and  profiles  are  compared  and  conclusions  drawn. 

The  governing  equations 

When  describing  multiphase  flow  phenomena  it  is  implied  that  more  than  one  phase  exist 
within  a  small  volume  at  any  particular  time.  This  view  rests  on  the  idea  of  time  and  space 
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averaging  [14,15]  and  is  equivalent  to  treating  each  phase  as  a  continuum  in  the  space  under 
consideration,  which  requires  choosing  a  proper  scale  with  regard  to  the  control  volume  used. 
For  a  multi-phase  system  the  equations  are  derived  over  a  representative  element  volume 
within  which  the  different  phases  are  present. 

Except  for  the  near  region  of  the  injector  where  the  spray  is  dense,  the  volume  fraction  of  the 
spray  is  low.  In  this  dilute  multi-phase  flow  regime,  interaction  between  droplets  can  be 
neglected.  Starting  from  the  Navier-Stokes  equations,  instantaneous  transport  equations  for 
the  gas  and  droplet  phase  can  be  derived  either  by  spatial,  temporal,  or  ensemble  averaging. 
However,  these  transport  equations  can  only  be  used  for  the  description  of  sprays  in  laminar 
gas  flows.  Since  combustors  generally  operate  in  the  turbulent  regime,  the  system  of 
equations  is  extended  by  introducing  turbulent  fluctuations  of  the  transport  quantities 
followed  by  Reynolds  averaging  of  the  equations.  For  the  gaseous  phase,  the  standard  k-e 
model  is  employed,  while  an  algebraic  model  based  on  a  Boussinesq  approach  approximates 
the  turbulence  terms  in  the  droplet  phase  transport  equations.  The  interacting  flow  fields  are 
described  by  the  transport  equations  presented  next. 


Gas  equations 

The  continuity,  momentum,  energy,  turbulence  kinetic  energy,  and  turbulence  dissipation  rate 
equations  for  the  gas  phase,  which  is  composed  of  two  species  namely  air  and  vapor,  in 
addition  to  the  mass  fraction  equation  of  the  fuel  vapor  in  the  gaseous  phase,  are  respectively 
written  as: 


^Kp,)+v-KugP,)=v- 


dt 

d_ 

dt 

d_ 

dt 


1,8 


Sc, 


Va, 


+  M 


vap,g 


(«.P«U.)  +  v-k  Wg)  =  ~^P  +  v-^  +  K  +F°  +  ^ud 

Kpa  ) + v-Kp*  v, ) = v-k  Nt,^k8 ) + k  -  pg£s ) + 5 
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dt 
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f  e  e 2A 
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dt 
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d_ 

dt 


vap,g 


)  +  V.(a::pi:u, 


)  “  ,eff^Yvap,g  )  +  M  mp  g 


(6) 


Droplet  balance  equations 

Droplet  evaporation  is  simulated  by  means  of  the  Uniform  Temperature  model  [16,17].  This 
computationally  effective  droplet  model  is  based  on  the  assumption  of  a  homogeneous 
internal  temperature  distribution  in  the  droplet  and  phase  equilibrium  conditions  at  the 
surface.  The  analytical  derivation  of  this  model  does  not  consider  contributions  to  heat  and 
mass  transport  through  forced  convection  by  the  gas  flow  around  the  droplet.  Forced 
convection  is  taken  into  account  by  means  of  two  empirical  correction  factors 

^correction  a,ul  KorrecUon)  [18]. 
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Correction  ,=  1  +  0.276  Re^Scf 


h  =  1 +0  276Re1/2  Prl/3 

correction  ,i  i  ^  VJ.Z.  /  U  n; 


Pt 


turb,d,i  P turb,g  , 

Pg  KS 

For  an  N-phase  flow,  the  volume  fractions  a(k)  are  characterized  by  the  condition: 


Pd,i  ^d,i 


(16) 

(17) 

(18) 


a 


w 


1 


(19) 


A--1 


The  ratio  of  the  turbulent  kinetic  energies  of  a  dispersed  (d)  and  gas  (g)  phase  is  calculated 
following  the  approach  in  [7,11]: 


1 


2_2 


1  +  fflf 


(20) 


where 


(21) 


kd  =  \^d-^d 

Since  in  general  the  droplets  do  not  follow  the  motion  of  the  surrounding  fluid  from  one  point 
to  another  it  is  expected  for  the  ratio  kd  /  kg  to  be  different  from  unity  and  varies  with  particle 

relaxation  time  t  and  local  turbulence  quantities.  Kramer  [11]  recommends  the  following 
equation  for  the  frequency  of  the  particle  response: 


1 


(O  =  — 
T 


fir 

V  3  s 


0/4 


T  = 


1  PdD 


V  *  J 
2 


(22) 


1 


0.687 


18  pg  v  1  +  0.133Re 

with  a  characteristic  macroscopic  length  scale  of  turbulence  given  by 
t  \3n 


(23) 
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Discretization  procedure 


A  review  of  the  above  differential  equations  reveals  that  they  are  similar  in  structure.  If  a 
typical  representative  variable  associated  with  phase  (k)  is  denoted  by  (f)(k> ,  the  general  fluidic 
differential  equation  may  be  written  as: 

d(aik)p(kYk)) 


dt 


■  + 


V.(a(i)p(t)u(iyt))=V.(a(i,r(i)V0(t))  +  a(i)0(t) 


(24) 


where  the  expression  for  Vkl  and  Q <k>  can  be  deduced  from  the  parent  equations. 

The  general  conservation  equation  (24)  is  integrated  over  a  finite  volume  to  yield: 

?( akk)  okk)(b(k)) 

- L_I_JrfQ  + JJv.(a<*yi,ul*>(i,)rfQ  =  jjV.(a(k)r(k)Vpm)dQ  +  \\awQmdQ. 

n  ri  n  n 


(25) 


Where  12  is  the  volume  of  the  control  cell.  Using  the  divergence  theorem  to  transform  the 
volume  integral  into  a  surface  integral,  replacing  the  surface  integrals  by  a  summation  of  the 
fluxes  over  the  sides  of  the  control  volume,  and  then  discretizing  these  fluxes  using  suitable 
interpolation  profiles  the  following  algebraic  equation  results: 

ApVpk)  =  +  Bkk)  (26) 

NB 

In  compact  form,  the  above  equation  can  be  written  as 

=  =  ^ - -p -  (27) 

Ap 

An  equation  similar  to  equation  (26)  is  obtained  at  each  grid  point  in  the  domain  and  the 
collection  of  these  equations  forms  a  system  that  is  solved  iteratively. 

The  discretization  procedure  for  the  momentum  equation  yields  an  algebraic  equation  of  the 
form: 

nf  =  HPp  [uw]  -  akk)D(pk)Vp  ( P )  (28) 

Furthermore,  the  phasic  mass-conservation  equation  can  be  viewed  as  a  phasic  volume 
fraction  equation,  which  can  be  written  as: 

a(k>  =  H  p[a(k>]  (29) 

or  as  a  phasic  continuity  equation  to  be  used  in  deriving  the  pressure  correction  equation: 


8t 


Q  +  Ap  [a(t)plt)u(t)JS]  =  a(k)M ik) 


(30) 
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where  the  A  operator  represents  the  following  operation: 

M®]=  X  ®/ 

f=nb(P ) 


(31) 


Pressure  correction  equation 

To  derive  the  pressure-correction  equation,  the  mass  conservation  equations  of  the  various 
fluids  are  added  to  yield  the  global  mass  conservation  equation  given  by: 


\Old 


St 


Q+Ap  (a(k)p<PkVk\s) 


=  0 


(32) 


Denoting  the  corrections  for  pressure,  density,  and  velocity  by  P',  n ,  and  p(k)  , 
respectively,  the  corrected  fields  are  written  as: 


P  =  P°  +  P',u(k)  =  u{kr  +  u  w,pw  =  p(k)°  +  p(ky 


(33) 


Combining  equations  (28),  (32),  and  (33),  the  final  form  of  the  pressure-correction  equation  is 
obtained  as  [19]: 


5^—  r{pk)°C^P;  + AP  [r(k)VarC(pk)p']-Ap  r(iyp(t)*(r(t)'D(t)VP')jS 

r«Vr-(r«>p<»)' 


~z  ■ 

k 


St 


n+A,[r(tVi)*tf(t)*] 


(34) 


The  corrections  are  then  applied  to  the  velocity,  density,  and  pressure  fields  using  the 
following  equations: 


Up )*  =  Up )o  -  ra)oDp  VpP',  P*  =  P°  +  P',  par  =  p(k)o  +  Cipk)P' 


(35) 


The  multiphase  flow  model 

In  the  multiphase  model,  displayed  schematically  in  Figure  1,  the  dispersed  phase  is 
subdivided  into  N  size  group  classes  where  each  class  is  treated  as  a  continuum  phase  in  the 
calculation.  Although  this  model  can  describe  a  real  multi-phase  flow  problem  however  the 
number  of  coupled  equations  associated  with  this  model  is  very  large.  The  number  of  size 
groups  should  therefore  be  limited  to  low  values  due  to  the  extensive  numerical  effort  of  the 
coupled  equations. 
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In  order  to  keep  track  of  the  droplet  size,  an  additional  diameter  equation  is  solved  for  every 
droplet  phase.  This  equation  is  derived  based  on  mass  conservation  and  is  given  by 


<Xi.iPi.id, )  +  V(<*i,Pi,iPdA )  =  V- 


%/,; 


V 


Pr, 


turb.d  ,1 


+  d;V. 


) 


f^turb,, 


d  ,i 


V  turb.d, i 


V  ad,i 


+  ~ dMdd  (36) 


This  multi-phase  model,  is  an  extension  of  the  work  performed  in  [20].  Results  generated 
using  this  model  are  the  baseline  against  which  results  generated  using  the  MUSIG  and  H- 
MUSIG  models  are  compared. 


Full  Multiphase  Solution  Procedure 

The  overall  solution  procedure  is  an  extension  of  the  single-phase  SIMPLE  algorithm  [21,22] 
into  multi-phase  flows  [19].  The  sequence  of  events  in  the  MCBA-SIMPLE  is  as  follows: 

1.  Solve  the  fluidic  momentum  equations  for  velocities. 

2.  Solve  the  pressure  correction  equation  based  on  global  mass  conservation. 

3.  Correct  velocities,  densities,  and  pressure. 

4.  Solve  the  fluidic  mass  conservation  equations  for  volume  fractions. 

5.  Solve  the  fluidic  scalar  equations  (k,  8,  T,  Y,  D,  etc. . .). 

6.  Return  to  the  first  step  and  repeat  until  convergence. 


The  Multi-Size-Group  model  (MUSIG) 

The  MUSIG  (MUltiple-SIze-Group)  model  [23,24]  represents  a  further  development  of  the 
multiphase  model  outlined  above.  The  MUSIG  model  was  originally  developed  for  the 
prediction  of  bubbles  in  water  and  has  never  been  used  for  the  prediction  of  mixing  and 
evaporation  of  liquid  droplets  in  a  stream  of  gas.  It  is  the  intention  of  this  work  to  extend  the 
applicability  of  MUSIG  to  such  configuration. 

As  a  first  step  for  decribing  the  MUSIG  model,  an  explanation  of  the  population  balance 
approach  is  provided.  This  is  followed  by  the  model  description,  the  equations  involved,  and 
the  break  up  and  coalescence  models  used. 
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Population  Balance  approach 

In  spray  modeling,  a  wide  range  of  particle  sizes  and  shapes  exist  at  every  point  in  the 
dispersed  phase,  which  makes  the  description  of  the  size  and  shape  of  the  droplets  very 
difficult.  This  difficulty  is  further  magnified  when  break  up  and  coalescence  occur  due  to 
their  great  influence  on  the  overall  performance  among  the  phases.  Therefore  it  is  essential  in 
modeling  spray  flows  to  use  a  formulation  that  take  into  account  the  different  size  distribution 
of  particles  in  addition  to  the  birth  and  death  processes  that  may  be  encountered  [25].  This  is 
accomplished  through  the  use  of  the  population  balance  approach. 

The  formulation  of  the  population  balance  could  be  either  based  on  the  Boltzmann  transport 
equation  or  on  continuum  mechanics  [26].  It  represents  the  transport  of  the  number  density  of 
the  fluid  through  the  space  taking  into  account  birth  and  death  of  particles  due  to  breakup  and 
coalescence.  The  fluid  particle  number  density  transport  equations  of  particles  having  volume 
Vi,  i.e.  group  size  i,  is  given  as  follows: 

^  +  V.(u «,)  =  S,+(S,1)[  (37) 

The  interaction  term  5,  represents  the  net  rate  of  change  in  the  number  density  distribution 
function,  due  to  particle  break-up  and  coalescence;  in  other  words  this  term  accounts  for 
both  the  birth  of  particles  of  size  /  (due  to  either  breakup  of  larger  particles  or  coalescence  of 
smaller  ones)  and  death  of  particles  of  size  i  (also  due  to  breakup  and  coalescence  but  this 
time  of  particles  of  size  i  which  goes  into  other  groups).  A  general  representation  of  these 
source  and  sink  terms  is  given  as 

Sj  =  Bm  -  DBi  +  BCi  -  Da  (38) 

Moreover  the  term  Sph  is  added  to  the  population  balance  equation  (Eq.  (37))  in  order  to 
account  for  phase  change  since  size  change  can  occur  because  of  birth  due  to  nucleation, 
condensation  and  evaporation  [27]. 

MUSIG  MODEL  DESCRIPTION 

In  this  model,  displayed  schematically  in  Figure  2,  the  dispersed  phase  is  decomposed  into  N 
size  groups  chosen  according  to  the  problem  at  hand  (say  N=10)  where  all  groups  are  treated 
as  one  continuum  phase  thus  moving  at  the  same  speed  and  N  continuity  or  population 
balance  equations  are  solved.  The  single  disperse  phase  is  thus  characterized  by  various  size 
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groups,  from  which  a  local  Sauter  mean  diameter  is  deduced.  An  accurate  determination  of 
the  droplet  Sauter  diameter  is  crucial  in  order  to  calculate  the  interfacial  area  based  on  which 
the  interphase  heat  and  mass  transfer  and  momentum  drags  could  be  evaluated. 

A  discretization  is  performed  in  order  to  evaluate  the  diameter  (d,)  of  each  size  group.  Any  of 
the  following  three  approaches  can  be  used  [24]:  (a)  the  equal  mass  discretization,  (b)  the 
equal  diameter  discretization,  (c)  The  geometric  mass  discretization.  In  all  three  methods  the 
minimum  and  maximum  diameters  of  the  polydispersed  fluid  are  set  before  beginning  the 
calculations. 

MUSIG  Model  formulation 

As  described  previously  a  continuity  equation  for  each  size  group  (a  population  balance 

equation)  is  solved,  but  it  is  assumed  that  all  droplets  move  with  the  same  velocity  so  that 

only  one  set  of  momentum  equations  for  the  dispersed  phase  has  to  be  solved. 

let  f  be  the  size  fraction  of  the  polydispersed  phase  which  appears  in  group  size  i,  pd  the 

density  of  the  dispersed  phase,  if  Equation  (38)  is  multiplied  by  m„  it  reduces  to  the 

following: 

m. 

p=  —  &  a  =  n  Q. 

d  Q  '  '  ' 

i 

ai  =  adf,  (39) 

Jj(pA/i)  +  v  '(PA/i11,)  =  s,  +(s  J, 

Since  this  model  assumes  that  all  particles  have  the  same  velocity,  «,  is  replaced  by  ud 
(velocity  of  the  dispersed  phase)  and  the  above  equation  reduces  to 

)  =  St+  (Sph );  (40) 

This  equation  has  the  form  of  the  transport  equation  of  a  scalar  variable  fi_  in  the  dispersed 
phase  in  which  the  source  term  S,  accounts  for  the  birth  of  droplets  of  size  i  due  to  breakup  of 
droplets  of  larger  size  and  coalescence  of  droplets  of  smaller  size  which  go  into  i  and  death  of 
droplets  of  size  i  due  to  both  break  up  and  coalescence  encountered  in  this  size  group. 
Therefore  the  sum  of  this  term  over  all  size  groups  is  equal  to  zero 

Xs.  =°  (41) 

i 
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If  summation  over  all  size  groups  is  performed,  an  overall  continuity  equation  for  the 
dispersed  phase  is  derived.  Performing  this  summation  and  noticing  that  ft  =  1  ,  the 

i=l,  AT 

overall  continuity  equation  for  the  dispersed  phase  is  found  as 
J;(PA/i)  +  =  Sl  -  fxMvaps 

^ -{Pd&dfl^d  )  =  ^2  _  flM  vap,g 

+  : 

;  (42) 


d_ 

dt 


( Pd^dfs  ) 


fNM 


vap.g 


Jt(P*ad  )  +  V-(PAUrf  )  =  v«P,S 

The  continuity  equation  is  therefore  used  to  solve  the  population  balance  equations  with  the 
source  term  to  be  connected  with  the  source  of  population  balance  equations.  It  is  important 
to  note  here  that  no  extra  scalar  equations  are  needed  for  the  population  balance  equation 
[28], 

After  solving  the  population  balance  equations,  the  droplet  sauter  diameter  which  represents 
an  average  representation  of  the  dispersed  phase  is  calculated,  and  the  dispersed  phase 
calculations  are  then  performed  assuming  the  monodispersed  phase  with  the  diameter  of  the 
phase  being  equal  to  the  local  droplet  Sauter  diameter  ds  based  on  the  values  of  the  size 
fractions  of  the  dispersed  phase/;  and  discrete  droplet  sizes  d,  given  by: 


As  in  the  multi-phase  model,  a  continuity  equation  for  each  size  group  is  solved,  but  it  is 
assumed  that  all  droplet  velocities  can  be  related  to  the  average  value  in  an  algebraic  manner, 
so  that  only  one  set  of  momentum  equations  for  the  gas  phase  is  solved.  This  single  disperse 
phase  is  characterized  by  various  size  groups,  from  which  a  local  Sauter  mean  diameter  is 
deduced.  An  accurate  determination  of  the  droplet  Sauter  diameter  is  crucial  in  order  to 
calculate  the  interfacial  area  based  on  which  the  interphase  heat  and  mass  transfer  and 
momentum  drags  could  be  evaluated. 
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The  MUSIG  model  essentially  reduces  the  multiphase  approach  described  above  back  to  a 
two-fluid  approach  with  one  velocity  field  for  the  continuous  phase  and  one  for  the  dispersed 
phase.  However,  the  continuity  equations  of  the  particle  size  groups  are  retained  and  solved  to 
represent  the  size  distribution.  With  this  approach,  it  is  possible  to  consider  a  larger  number  of 
particle  size  groups  (say  10,  20  or  even  30  particle  phases)  to  give  a  better  representation  of 
the  size  distribution. 


Break  up  and  coalescence  models 

Various  models  have  been  implemented  for  the  calculation  of  the  break-up  and  coalescence 
rates.  The  overall  form  of  these  models  is  similar,  with  different  correlations  used  for 
coalescence  and  break-up  depending  on  the  nature  of  the  flow  [27] . 

The  specific  models  that  has  been  implemented  for  the  break-up  and  coalescence  rates  are  due 
to  Luo  and  Svendsen  [12]  and  Coulaloglou  and  Tavlarides  [13]. 


Luo  and  Svendsen  Breakup  Model 


The  net  source  to  group  i  due  to  breakup  (using  any  break  up  model)  is: 

f  \ 


A  -  =  />.«  : 


I"  -'.-/I". 

V  i>'  j<‘  J 


(44) 


The  break-up  rate  BtJ  is  assumed  to  be  a  function  of  the  break-up  fraction  B\j  as  follows: 


B,  =  B',  j  df„  (45) 

fav 

This  model  is  based  on  the  following  five  main  assumptions: 

1-  The  turbulence  is  assumed  to  be  isotropic. 

2-  Binary  breakage  of  fluid  particles  is  only  considered  (turbulent  breakage  and  Shear 
breakage).  The  breakage  volume  fraction  is  assumed  to  be  given  by 


_mL_d1_ 

JBV  ,4 


m . 


d 


(46) 


Where  dj  is  the  diameter  of  a  mother  particle  splitting  into  two  particles  of  sizes  d,  and 
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3-  The  occurrence  of  breakup  is  determined  by  the  energy  level  of  the  arriving  eddy.  The 
particle  oscillation  frequency  is  larger  than  the  arrival  frequency  of  eddies.  This 
implies  that  eddies  affect  the  particles  independently  in  a  way  that  once  an  eddy  of 
sufficiently  high  energy  arrives,  the  particle  will  break. 

4-  Eddies  of  size  scale  smaller  than  or  equal  to  the  particle  diameter  induce  particle 
oscillation  and  dominate  the  deformation  of  the  particles  in  the  flow  field  while  larger 
eddies  are  responsible  for  the  translatory  motion  of  particles. 

The  break-up  frequency  of  a  mother  particle  with  size  dj  splitting  into  two  particles  of  sizes  d, 

and  (dj  -  d-  )  is  given  by: 


i 


where  Fb  is  added  as  a  calibration  factor  of  the  model,  sg  is  the  continuous  phase  eddy 
dissipation  energy,  a  is  the  surface  tension,  /?=2,  and  <f  is  the  dimensionless  size  of  eddies  in 
the  inertial  subrange  of  isotropic  turbulence.  The  lower  limit  of  the  integration  is  given  by 


£  . 
^min 


=  11.4 


a 

d. 

J 


and  the  Kolmogorov  micro-scale  q  is  given  by: 


n  = 


,1/4 


(48) 


(49) 


COULALOGLOU  AND  TAVLARIDES  COALESCENCE  MODEL 
Collisions  may  happen  due  to  various  mechanisms  however  this  model  only  considers  the 
collisions  of  droplets  due  to  turbulence,  buoyancy  and  laminar  shear. 

The  net  source  to  group  i  due  to  coalescence  is  given  as 


(l 

Ba~Da  =(Aa4  rXSCj//, 


lit  +  mk 
">i»h 


Xj«  ~ 


1^: 


mjj 


(50) 
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where  Q  is  the  specific  coalescence  rate  between  groups  i  and  j  and  Xjki  represents  the 
fraction  of  mass  due  to  coalescence  between  groups  j  and  k  which  goes  into  i  and  is  written  as 


X 


jki 


(mj  +  mk 
mi  —  mi  j 
m,+i  -  (m j  +  mk ,) 

mi+i  ~  m, 

1 

0 


if  mi_l<(mj  +  mk)<mi 

if  mi<(mj+mk)<mi+l 

if  nij  +  mk  >  rnmax  =  m, 
otherwise 


'LX  lk<  =  ]  for  all  j ,  k 


(51) 


(52) 


When  summed  over  all  size  groups,  the  net  source  due  to  coalescence  is  zero. 

The  coalescence  rate  Q  of  the  dispersed  phase  drops  in  a  turbulent  flow  field  is  described  as 
the  product  of  the  collision  frequency  Q)c(dj,dj}  and  the  corresponding  coalescence 
efficiency  A  (//,,^()  as  is  written  as 


C,  =  a 


(53) 


In  the  Coulaloglou  and  Tavlarides  [13]  model  it  is  assumed  that  the  mechanism  of  collision  in 
a  locally  isotropic  flow  field  is  analogous  to  collisions  between  molecules  as  in  kinetic  theory 
of  gases,  the  collision  frequency  between  two  drops  with  volumes  Q,  and  Q,  are  expressed  as 


1/3 

‘°M-dj)  =  Kcffx-(di  +djf  Uf+dff 


(54) 


The  coalescence  efficiency  is  based  on  the  film  drainage  mechanism  where  the  drops  are 
considered  to  cohere  together  and  be  prevented  by  coalescence  by  a  film  of  continuous  phase 
trapped  between  them.  The  coalescence  efficiency  as  suggested  by  Coulaloglou  and 
Tavlarides  is  given  as 


A(CW,): 


exp 
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VcP,£c 
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d ,  +  d 
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(55) 
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MUISG  Solution  Procedure 

The  overall  two-phase  MUSIG  solution  procedure  is  as  follows: 

1.  Solve  the  fluidic  momentum  equations  for  the  gas  and  droplet  phase  for  velocities. 

2.  Solve  the  pressure  correction  equation  based  on  global  mass  conservation. 

3.  Correct  velocities,  densities,  and  pressure. 

4.  Solve  the  fluidic  mass  conservation  equations  for  droplet  and  air  volume  fractions. 

5.  Calculate  sources  due  breakup  and  coalescence. 

6.  Solve  the  population  fraction  equations. 

7.  Calculate  Sauter  diameter. 

8.  Solve  the  fluidic  scalar  equations  (k,  8,  T,  Y). 

9.  Return  to  the  first  step  and  repeat  until  convergence. 

The  Heterogeneous  MUSIG  Model  (H-MUSIG) 

With  the  MUSIG  model,  it  is  possible  to  consider  a  larger  number  of  particle  size  groups  to 
give  a  better  representation  of  the  size  distribution.  The  shortcoming  of  this  approach 
however  [29,  30],  is  related  to  the  droplet  groups  common  velocity.  It  is  well  known  that 
larger  droplets  do  not  follow  the  flow  and  smaller  droplets  do.  By  considering  one  average 
velocity  for  the  droplets,  a  stratification  of  droplet  sizes  from  normal  fuel  injection  occurs 
with  larger  droplets  penetrating  further  into  the  flow.  The  larger  droplets  transport  more  fuel 
mass  than  may  be  expected.  To  alleviate  this  problem,  it  is  proposed  to  extend  the 
(homogeneous)  MUSIG  model  into  a  Multi-phase  MUSIG  model  or  a  Heterogeneous  MUSIG 
model  (H-MUSIG).  In  the  extended  model,  rather  than  assigning  one  velocity  for  all  droplet 
groups,  classes  of  groups  will  be  considered  with  droplet  groups  in  a  class  sharing  the  same 
velocity.  This  suggested  approach  ,  displayed  schematically  in  Figure  3,  could  be  seen  as  a 
blend  between  a  full  multi-phase  approach  and  a  two-phase  approach.  If  a  group  is  composed 
of  one  droplet  class,  then  the  full  multi-phase  approach  is  obtained,  whereas  if  a  group  is 
composed  of  all  droplet  sizes,  then  the  original  MUSIG  is  recovered. 
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H-MUSIG  Model  Description 


In  this  model,  the  dispersed  phase  is  divided  into  N  fields,  each  allowing  an  arbitrary  number 
of  sub-size  classes  [29,30].  Therefore  N  velocity  fields  are  to  be  solved  where  each  field  is 
subdivided  into  Mi  size  groups  moving  at  the  same  mean  algebraic  velocity.  Mi  population 
fraction  equations  need  to  be  solved  along  with  one  set  of  momentum  equations  for  each 
field. 

The  sub-size  groups  are  discretized  in  order  to  evaluate  the  diameter  (d,)  of  each  size  group. 
The  droplet  average  diameter  for  each  sub-size  class  is  calculated  using  one  of  the  group  size 
algorithms  described  above.  The  lower  and  upper  boundaries  for  each  field  should  be  set  in 
order  for  the  average  diameter  to  be  calculated. 


H-MUSIG  Model  Formulation 

As  described  earlier,  a  continuity  or  a  population  balance  equation  for  each  sub- size  group  is 
solved.  Let  /,  /e[l,N]  be  a  group  field  number;  the  volume  fractions  of  sub-classes,  groups  and 
the  dispersed  phase  are  related  as  follows: 

~  &dfi  ~  ^F,  f Fj  ,i  (56) 

where/  is  the  size  fraction  of  the  polydispersed  phase  which  appears  in  group  sub-size  i,  fij  is 
the  population  fraction  of  sub-group  size  i  in  the  field  /,  pd  the  density  of  the  dispersed  phase, 
ai  is  the  volume  fraction  of  sub-group  size  i  and  aF  is  the  volume  fraction  of  field  /.  If  Eq. 

(40)  is  multiplied  by  m„  the  population  fraction  equation  for  each  sub-size  group  in  class  field 
l  reduces  to 


dt 


[Pd^F,  f h),!  )  “I"  /)-  )  —  S Ft ,,  +  f ph 

OCj 


(57) 


Summing  over  all  sub-classes  in  field  l  and  applying  the  following  additional  relations 


CCr, 


i=i  ad 


a E 


a. 


(58) 


The  continuity  equation  for  each  field  or  the  volume  fraction  equation  is  found  as 


dt 


{PdaF, )  +  V .[pdaFuFi )  —  +  S ph 

OCj 


(59) 
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Additional  relations  and  constraints  are  further  added  to  this  model: 

N  NxM  mf,  NxM 

ad  =  IX  =  X  ai  IX  = 1  X  fi  = 1  (60) 

/  =  1  !  =  1  1  =  1  (  =  1 

The  source  SF(  is  the  net  rate  at  which  mass  accumulates  in  sub-size  group  i  due  to 

coalescence  of  particles  of  smaller  size  than  i  in  all  field  and  breakup  of  particles  of  larger 
sizes  than  i  in  all  fields.  Mathematicall  SF  is  given  by 

SF,  ,i  =  Bi,B  -  Di.B  +  Bi,c  ~  D,  c  (6 1 ) 

Here  the  birth  and  death  rates  are  treated  using  specific  breakup  and  coalescence  models  since 
they  represents  rates  in  the  dispersed  phase  for  all  sub-classes  with 


mf,  N 

S,,=2X,  &  2>„=°  (62) 

i=l  1=1 

The  overall  continuity  of  the  dispersed  phase  is  derived  by  summing  over  all  groups  using  the 
above  mentioned  relation  and  is  found  to  be 


d_ 

dt 


(ft«J  +  V.(pA  uF; 


(63) 


After  solving  the  population  balance  equations  for  each  sub-size  group,  the  droplet  sauter 
diameter  that  represents  an  average  representation  of  the  dispersed  phase  field  l  is  calculated, 
and  the  phase  field  calculations  are  then  performed  assuming  the  monodispersed  phase  for 
each  field  with  the  diameter  of  the  phase  field  being  equal  to 


d 


s,F, 


1 


MF/ 

1 


(64) 


For  each  field,  the  momentum  equation  that  need  to  be  solved  (since  the  velocities  differ)  is 
given  by 

|"  (  Pd  aF,  UF,  )  +  V‘  (  Pd  aF,  UF,  UF,  )  =  -aF,VP  +  (65) 

V.QJf|  /J-Iurbj  (Vn Fi  +  j  +  cxFi  pdg  +  FFi  +  M  Fi  +  SM^ 

where  M  F  is  the  source  term  associated  with  momentum  due  to  phase  change,  is  the 
momentum  transfer  between  the  dispersed  field  /  phase  and  the  continuous  phase,  and  is 
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the  source  term  due  to  the  transfer  of  momentum  between  different  velocity  groups  due  to 
breakup  and  coalescence  processes  leading  to  the  formation  of  particles  belonging  to  other 
groups. 

H-MUSIG  Model  Solution  Strategy 

The  solution  algorithm  of  the  system  is  as  follows: 

1.  The  momentum  equations  for  all  phases  are  solved  in  order  to  calculate  the  velocities 
of  the  different  phase  fields  and  the  continuous  phase. 

2.  The  Global  continuity  equation  derived  by  summing  continuity  for  all  phases  is  than 
solved  to  calculate  the  pressure,  which  is  the  same  for  all  phases. 

3.  The  continuity  equations  or  volume  fraction  equations  are  then  solved  to  calculate  the 
volume  fraction  of  the  dispersed  phases  (note  that  the  volume  fraction  equations  of  the 
various  phase  fields  include  a  source  term  due  to  particles  break  up  and  coalescence). 

4.  Calculate  sources  due  breakup  and  coalescence. 

5.  Solve  the  population  fraction  equations  for  the  groups  of  each  droplet  phase. 

6.  Calculate  Sauter  diameters  of  the  various  phases. 

7.  Solve  the  fluidic  scalar  equations  for  all  phases  (k,  e,  T,  Y). 

8.  Return  to  the  first  step  and  repeat  until  convergence. 
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Comparison  Between  MUSIC  and  H-MUSIG 
A  summary  of  the  differences  between  MUSIG  and  H-MUSIG  is  shown  in  the  table  below. 


Table  1  Comparison  between  homogeneous  and  inhomogeneous  model 


H-MUSIG 

MUSIG 

Dispersed  phase  accounts  for  N  phase 

fields  or  classes 

Dispersed  phase  accounts  for  one  phase 

field  or  class 

Each  class  is  subdivided  into  M  groups 

The  class  is  divided  into  N  groups 

N  sets  of  momentum  equations  needs  to 

be  solved 

1  set  of  momentum  equations  needs  to  be 

solved 

Source  term  due  to  break  up  and  coalescence  resulting  in  the  birth  or  death  of  particle  of 

size  i  in  the  population  balance  equations  is  the  same 

The  volume  fraction  equation  for  each 

class  1  includes  a  source  term  associated 

with  the  particle  break  up  and  coalescence 

contributions  (termSF  ) 

No  Contributions  for  particle  breakup  and 

coalescence  in  the  volume  fraction 

equation  since  no  other  dispersed  class  is 

present 

Interphase  momentum  and  energy  transfer 

term  due  to  the  presence  of  various 

velocity  and  energy  fields  are  added  to  the 

momentum  and  energy  equations, 

respectively. 

Results  and  Discussion 


In  what  follows,  solutions  to  three  two-dimensional  multi-phase  flow  problems  are  presented 
and  discussed.  The  physical  situations  for  these  problems  are  displayed  in  Figure  4.  Figure 
4(a)  represents  a  rectangular  duct  in  which  air  enters  with  a  uniform  free  stream  velocity  U, 
while  fuel  (Kerosene  is  used  in  all  computations  due  to  the  unavailability  of  the  physical 
properties  of  JP-7  fuel  to  the  author)  mixed  with  air  is  injected  through  a  nozzle  2  mm  in 
diameter  in  the  stream- wise  direction  through  120°  angle.  Figure  4(b)  is  similar  to  Figure  4(a) 


26 


with  the  exception  of  the  domain  being  cylindrical  (axi- symmetric)  and  as  such  only  half  the 
cylindrical  enclosure  is  considered.  Figure  4(c)  represents  the  same  rectangular  duct  displayed 
in  Figure  4(a)  with  fuel  being  injected  in  the  cross-stream  direction.  The  length  of  the  domain 
is  L  and  its  width  is  W  (W=L/4).  Turbulent  flow  results  are  generated  using  the  k-e 
turbulence  models.  Solutions  for  the  above  configurations  are  generated  using  the  various 
methodology  and  results  are  compared. 

Streamwise  Injection  in  a  Rectangular  Domain 

The  physical  domain  depicted  in  Fig.  4(a)  is  subdivided  into  120x102  non-uniform  control 
volumes.  The  length  L  of  the  domain  is  1  m.  The  fuel  is  injected  with  a  velocity  of  magnitude 
30  m/s  through  12  uniform  control  volumes  (each  of  width  .001/12  m)  at  different  injection 
angles  (varying  uniformly  from  -60°  to  60°  as  shown  in  Fig.  4(a)).  In  order  to  show  the 
applicability  of  the  solutions  procedures  for  fluid  flowing  at  all  speeds,  results  for  this 
configuration  are  generated  for  fuel  injection  in  both  subsonic  and  supersonic  streams. 

Evaporation  and  mixing  of  fuel  droplets  in  a  subsonic  stream 

For  the  physical  situation  depicted  in  Fig.  4(a),  the  Mach  number  and  temperature  of  the  air  at 
inlet  to  the  domain  are  taken  to  be  0.2  (Mairinlet=0.2)  and  700  K,  respectively.  The  mixture  of 
air  and  droplets  are  injected  into  the  domain  with  a  temperature  of  350  K  with  the  volume 
fraction  of  Kerosene  in  the  injected  air-fuel  mixture  being  0.1.  The  velocity  of  the  injected 
mixture  is  set  at  30  m/s  with  the  angle  of  injection  varying  from  -60°  to  60°.  With  this 
velocity  profile  and  volume  fraction  a  total  of  1.8327  Kg/s/m  of  fuel  are  injected  into  the 
domain. 

Full  multiphase  results  are  generated  using  5  droplet  phases  of  sizes  of  60  (tm,  80  pm,  100  pm, 
120  pm,  and  140  pm  with  their  inlet  volume  fractions  being  0.0125,  0.0225,  0.03,  0.0225,  and 
0.0125  respectively.  For  the  MUSIG  model,  the  droplet  phase  is  divided  into  10  groups  using 
the  equal  diameter  discretization  with  the  diameter  of  the  smallest  droplet  set  to  55  pm  and  the 
increment  to  10  pm  with  population  fractions  of  0.05,  0.075,  0.1,  0.125,  0.15,  0.15,  0.125, 
0.1,  0.075,  and  0.05,  respectively.  For  the  H-MUSIG  model,  two  droplet  phases  are 
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considered,  each  divided  into  five  group  sizes  discretized  using  the  equal  diameter 
discretization.  The  diameters  and  population  fractions  of  the  various  groups  are  similar  to 
those  used  with  MUSIG.  An  important  parameter  for  comparison  is  the  the  percent  of  the 
injected  fuel  that  has  evaporated  into  the  gas  field.  These  percentages  are  found  to  be  30.83%, 
27.78%,  and  27.27%  for  the  full  multiphase  method,  the  H-MUSIG  model,  and  the  MUSIG 
model  respectively. 

Comparison  of  results  obtained  using  the  various  techniques  are  presented  in  Figures  5 
through  10.  Figure  5  displays  the  velocity  fields  for  some  of  the  droplet  phases.  Figure  5(a)- 
5(c)  depicts  the  velocity  vectors  for  droplet  phases  1  (60  (im  in  diameter),  3  (100  pm  in 
diameter),  and  5  (140  pm  in  diameter)  using  the  full  multiphase  approach.  H-MUSIG  droplet 
velocity  vectors  are  presented  in  Figures  5(d)  and  5(e),  while  Figure  5(f)  shows  the  droplet 
vector  field  predicted  using  MUSIG.  Both  the  multiphase  and  the  H-MUSIG  results  reveal  a 
larger  droplet  peneteration  with  increasing  droplet  diameter,  which  is  physically  correct  as 
larger  particles  possess  higher  inertia  and  are  more  capable  of  penetrating  into  the  domain  as 
compared  to  smaller  ones,  which  align  faster  with  the  flow  field.  Further  the  path  of  the 
droplets  predicted  by  MUSIG  is  between  the  trajectories  of  the  smaller  and  larger  droplet 
phases  predicted  by  H-MUSIG  (Compare  Fig.  5(f)  against  Figures  5(d)  and  5(e)).  The  same  is 
true  for  H-MUSIG  and  the  full  multiphase  results  (compare  Figure  5(d)  against  Figures  5(a) 
and  5(b);  and  Figure  5(e)  against  Figures  5(b)  and  5(c)). 

Comparisons  in  the  form  of  contour  maps  of  the  gas  phase  volume  fraction,  vapor  mass 
fraction,  pressure  field,  and  gas  temperature  fields  generated  by  the  various  methods  are 
depicted  in  Figures  6,  7,  8,  and  9  respectively.  As  can  be  seen,  the  overall  structure  of  the 
various  fields  are  similar  even  though  there  exist  some  variations  in  the  details.  The  volume 
fraction  fields  (Figure  6)  reflect  the  droplet  velocity  fields  displayed  in  Figure  5,  with  the 
droplet  volume  fraction  decreasing  as  droplets  move  in  the  domain.  As  expected,  the  vapor 
mass  fraction  (Figure  7)  in  the  gas  phase  maximizes  at  exit  from  the  domain  with  the  full 
multiphase  results  showing  less  evaporation  in  the  region  around  the  centerline  of  the  domain. 
Pressure  and  temperature  field  maps  presented  in  Figures  8  and  9  show  similar  profiles  with 
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slight  variation  in  levels.  This  is  further  revealed  in  the  comparison  of  profiles  presented  in 
Figure  10.  In  this  figure,  u- velocity,  gas  temperature,  pressure,  and  vapor  mass  fraction 
profiles  across  the  domain  at  x=0.5  generated  using  the  various  algorithms  are  compared.  The 
gas  u-velocity  profiles  (Figure  10(a))  obtained  by  the  various  methods  are  nearly  coincident. 
The  temperature  profiles  however  (Figure  10(b)),  show  some  differences  in  the  region  around 
the  centerline  of  the  domain,  with  H-MUSIG  predicting  the  lowest  gas  temperature.  This  is  in 
line  with  the  vapor  mass  fraction  contours  presented  earlier  and  the  vapor  mass  fraction 
profiles  in  Figure  10(d),  which  reveal  higher  vapor  mass  fraction  around  the  centerline 
predicted  by  H-MUSIG  and  consequently  lower  gas  temperature.  Moreover,  profiles  in  Figure 
10(d)  show  that  full  multiphase  results  are  close  to  H-MUSIG  results  in  area  away  from  the 
centerile  and  are  close  to  MUSIG  results  in  the  region  around  the  centerline.  Pressure  profiles 
presented  in  figure  10(c)  indicate  lower  level  values  with  the  full  multiphase  approach. 
Values  obtained  by  MUSIG  and  H-MUSIG  are  very  close. 

Evaporation  and  mixing  of  fuel  droplets  in  a  supersonic  stream 

For  the  physical  situation  depicted  in  Fig.  4(a),  the  Mach  number  and  temperature  of  the  air  at 
inlet  to  the  domain  are  taken  to  be  2  (Mair  inlet=2)  and  700  K,  respectively.  The  mixture  of  air 
and  droplets  are  injected  into  the  domain  with  a  temperature  of  350  K  with  the  volume 
fraction  of  Kerosene  in  the  injected  air-fuel  mixture  being  0.1.  The  velocity  of  the  injected 
mixture  is  set  at  30  m/s  with  the  angle  of  injection  varying  from  -60°  to  60°.  With  this 
velocity  profile  and  volume  fraction  a  total  of  1.8327  Kg/s/m  of  fuel  are  injected  into  the 
domain. 

Full  multiphase  results  are  generated  using  5  droplet  phases  of  sizes  of  60  (tm,  80  pm,  100  pm, 
120  pm,  and  140  pm  with  their  inlet  volume  fractions  being  0.0125,  0.0225,  0.03,  0.0225,  and 
0.0125  respectively.  For  the  MUSIG  model,  the  droplet  phase  is  divided  into  10  groups  using 
the  equal  diameter  discretization  with  the  diameter  of  the  smallest  droplet  set  to  55  pm  and  the 
increment  to  10  pm  with  population  fractions  of  0.05,  0.075,  0.1,  0.125,  0.15,  0.15,  0.125, 
0.1,  0.075,  and  0.05,  respectively.  For  the  H-MUSIG  model,  two  droplet  phases  are 
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considered,  each  divided  into  five  group  sizes  discretized  using  the  equal  diameter 
discretization.  The  diameters  and  population  fractions  of  the  various  groups  are  similar  to 
those  used  with  MUSIG.  An  important  parameter  for  comparison  is  the  the  percent  of  the 
injected  fuel  that  has  evaporated  in  the  domain.  These  percentages  are  found  to  be  9.69%, 
12%,  and  12.38%  for  the  full  multiphase  method,  the  H-MUSIG  model,  and  the  MUSIG 
model  respectively. 

Comparison  of  results  obtained  using  the  various  techniques  are  presented  in  Figures  11 
through  16.  Figure  11  displays  the  velocity  fields  for  some  of  the  droplet  phases.  Figure 
1 1(a)- 1 1(c)  depicts  the  velocity  vectors  for  droplet  phases  1  (60  pm  in  diameter),  3  (100  pm 
in  diameter),  and  5  (140  pm  in  diameter)  using  the  full  multiphase  approach.  H-MUSIG 
droplet  velocity  vectors  are  presented  in  Figures  11(d)  and  11(e),  while  Figure  11(f)  shows 
the  droplet  vector  field  predicted  using  MUSIG.  Due  to  the  small  velocity  by  which  the  fuel  is 
injected  (30  m/s)  as  compared  to  the  gas  velocity  (1061  m/s)  the  spread  is  much  lower  than 
the  subsonic  case.  Nevertheless,  both  the  multiphase  and  the  H-MUSIG  results  reveal  a  larger 
droplet  peneteration  with  increasing  droplet  diameter,  which  is  physically  correct  as  larger 
particles  possess  higher  inertia  and  are  more  capable  of  penetrating  into  the  domain  as 
compared  to  smaller  ones,  which  align  faster  with  the  flow  field.  Further  the  path  of  the 
droplets  predicted  by  MUSIG  is  between  the  trajectories  of  the  smaller  and  larger  droplet 
phases  predicted  by  H-MUSIG  (Compare  Fig.  11(f)  against  Figures  11(d)  and  11(e)).  The 
same  is  true  for  H-MUSIG  and  the  full  multiphase  results  (compare  Figure  11(d)  against 
Figures  11(a)  and  11(b);  and  Figure  11(e)  against  Figures  11(b)  and  11(c)). 

Comparisons  in  the  form  of  contour  maps  of  the  gas  phase  volume  fraction,  vapor  mass 
fraction,  pressure  field,  and  gas  temperature  fields  generated  by  the  various  methods  are 
depicted  in  Figures  12,  13,  14,  and  15  respectively.  As  can  be  seen,  the  overall  structure  of  the 
various  fields  are  similar  even  though  there  exist  some  variations  in  the  details.  The  volume 
fraction  fields  (Figure  12)  reflect  the  droplet  velocity  fields  displayed  in  Figure  11,  with  the 
droplet  volume  fraction  decreasing  as  droplets  move  in  the  domain.  As  expected,  the  vapor 
mass  fraction  (Figure  13)  in  the  gas  phase  maximizes  at  exit  from  the  domain  with  the  full 
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multiphase  results  showing  less  evaporation  in  the  region  around  the  centerline  of  the  domain. 
Pressure  and  temperature  field  maps  presented  in  Figures  14  and  15  show  similar  profiles 
with  slight  variation  in  levels.  This  is  further  revealed  in  the  comparison  of  profiles  presented 
in  Figure  16.  In  this  figure,  u- velocity,  gas  temperature,  pressure,  and  vapor  mass  fraction 
profiles  across  the  domain  at  x=0.5  generated  using  the  various  algorithms  are  compared.  The 
gas  u-velocity  profiles  (Figure  16(a))  obtained  by  the  various  methods  are  nearly  coincident. 
The  temperature  profiles  however  (Figure  16(b)),  show  a  slight  difference  in  the  region 
around  the  centerline  of  the  domain,  with  profiles  predited  by  MUSIG  and  H-MUSIG  being 
almost  coincident.  Pressure  profiles  presented  in  figure  16(c)  show  similar  variations  with 
values  obtained  with  the  full  multiphase  method  being  slightly  higher.  Moreover,  profiles  in 
Figure  16(d)  show  that  H-MUSIG  results  are  closer  to  the  full  multiphase  results  with  the 
differences  between  the  various  profiles  being  in  the  region  around  the  centerline  and  the 
largest  evaporation  rate  being  predicted  by  the  full  multiphase  method. 

Cross-Stream  Injection  Into  a  Supersonic  Stream  in  a  Rectangular  Domain 

The  physical  domain  depicted  in  Fig.  4(b)  is  subdivided  into  130x70  non-uniform  control 
volumes.  The  length  L  of  the  domain  is  1.1  m.  The  fuel  is  injected  through  12  uniform  control 
volumes  (each  of  width  .001/12  m)  from  two  nozzles  located  on  the  lower  and  upper  walls  at 
10  cm  from  the  inlet.  The  Mach  number  and  temperature  of  the  air  at  inlet  to  the  domain  are 
taken  to  be  2  (Mair  inlet=2)  and  700  K,  respectively.  The  mixture  of  air  and  droplets  are  injected 
into  the  domain  at  a  temperature  of  350  K  with  the  volume  fraction  of  Kerosene  in  the 
injected  air-fuel  mixture  being  0.01.  The  velocity  of  the  injected  mixture  is  150  m/s  and  the 
angle  of  injection  is  60°.  With  this  velocity  and  volume  fraction  a  total  of  2.34  Kg/s/m  of  fuel 
are  injected  into  the  domain. 

Full  multiphase  results  are  generated  using  5  droplet  phases  of  sizes  of  60  (tm,  80  pm,  100  pm, 
120  pm,  and  140  pm  with  their  inlet  volume  fractions  being  0.0025,  0.0045,  0.006,  0.0045, 
and  0.0025,  respectively.  For  the  MUSIG  model,  the  droplet  phase  is  divided  into  10  groups 
using  the  equal  diameter  discretization  with  the  diameter  of  the  smallest  droplet  set  to  55  pm 
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and  the  increment  to  10  (tm  with  population  fractions  of  0.05,  0.075,  0.1,  0.125,  0.15,  0.15, 
0.125,  0.1,  0.075,  and  0.05,  respectively.  For  the  H-MUSIG  model,  two  droplet  phases  are 
considered,  each  divided  into  five  group  sizes  discretized  using  the  equal  diameter 
discretization.  The  diameters  and  population  fractions  of  the  various  groups  are  similar  to 
those  used  with  MUSIG.  The  percentages  of  the  injected  fuel  that  has  evaporated  into  the  gas 
field  as  predicted  by  the  various  methods  are  found  to  be  10.64%,  13.25%,  and  13.41%  for 
the  full  multiphase  method,  the  H-MUSIG  model,  and  the  MUSIG  model  respectively. 
Comparison  of  results  obtained  using  the  various  techniques  are  presented  in  Figures  17 
through  22.  Figure  17  displays  the  velocity  fields  for  some  of  the  droplet  phases.  Figure 
17(a)-  17(c)  depicts  the  velocity  vectors  for  droplet  phases  1  (60  (tm  in  diameter),  3  (100  (tm 
in  diameter),  and  5  (140  (im  in  diameter)  using  the  full  multiphase  approach.  H-MUSIG 
droplet  velocity  vectors  are  presented  in  Figures  17(d)  and  17(e),  while  Figure  17(f)  shows 
the  droplet  vector  field  predicted  using  MUSIG.  As  in  the  streamwise  injection  case,  the 
multiphase  and  the  H-MUSIG  results  reveal  a  relatively  larger  droplet  peneteration  with 
increasing  droplet  diameter.  Further  the  path  of  the  droplets  predicted  by  MUSIG  is  between 
the  trajectories  of  the  smaller  and  larger  droplet  phases  predicted  by  H-MUSIG  (Compare 
Fig.  17(f)  against  Figures  17(d)  and  17(e)).  The  same  is  true  for  H-MUSIG  and  the  full 
multiphase  results  (compare  Figure  17(d)  against  Figures  17(a)  and  17(b);  and  Figure  17(e) 
against  Figures  17(b)  and  17(c)). 

Comparisons  in  the  form  of  contour  maps  of  the  gas  phase  volume  fraction,  vapor  mass 
fraction,  pressure  field,  and  gas  temperature  fields  generated  by  the  various  methods  are 
depicted  in  Figures  18,  19,  20,  and  21  respectively.  The  volume  fraction  fields  (Figure  18) 
reflect  the  droplet  velocity  fields  displayed  in  Figure  17,  with  the  droplet  volume  fraction 
decreasing  as  droplets  move  in  the  domain.  As  expected,  the  vapor  mass  fraction  (Figure  19) 
in  the  gas  phase  maximizes  at  exit  from  the  domain.  With  the  injection  velocity  used,  droplets 
are  not  capable  to  penetrate  into  the  central  core  of  the  domain.  Pressure  and  temperature  field 
maps  presented  in  Figures  20  and  21  show  similar  profiles  with  slight  variation  in  levels. 
Interaction  of  both  injected  streams  via  pressure  is  obvious. 
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In  Figure  22,  the  u-velocity,  gas  temperature,  pressure,  and  vapor  mass  fraction  profiles 
across  the  domain  at  x=0.6  (i.e.  at  0.5  m  from  the  nozzles)  generated  using  the  various 
algorithms  are  compared.  As  in  the  previous  cases,  the  gas  u-velocity  (Figure  22(a)) 
component  and  gas  temperature  (Figure  22(b))  profiles  are  nearly  coincident.  Moreover, 
pressure  valued  (Figure  22(c))  are  closer  then  in  the  previous  cases  with  MUSIG  and  H- 
MUSIG  profiles  being  on  top  of  each  others.  Furthermore,  vapor  mass  fraction  profiles 
displayed  in  Figure  22(d)  indicate  that  H-MUSIG  results  are  closer  to  the  full  multiphase 
results  with  maximum  differences  between  the  various  profiles  occuring  in  the  region  close  to 
the  lower  and  upper  walls. 

Streamwise  Injection  Into  a  Supersonic  Stream  in  a  Cylindrical  Domain 

The  physical  domain  depicted  in  Figure  4(c)  representing  the  front  view  of  the  upper  half  of  a 
cylindrical  duct  of  length  1  m,  is  subdivided  into  120x60  non-uniform  control  volumes.  The 
fuel  is  injected  through  12  uniform  control  volumes  (each  of  width  .001/12  m)  from  a  nozzle 
located  at  the  center  of  the  domain,  as  shown  in  Figure  4(c).  By  neglecting  body  forces,  the 
problem  can  be  solved  as  an  axi-symmetric  two-dimensional  one.  The  Mach  number  and 
temperature  of  the  air  at  inlet  to  the  domain  are  taken  to  be  2  (Mair  inlet=2)  and  700  K, 
respectively.  The  mixture  of  air  and  droplets  are  injected  into  the  domain  at  a  temperature  of 
350  K  with  the  volume  fraction  of  Kerosene  in  the  injected  air-fuel  mixture  being  0.5.  The 
velocity  of  the  injected  mixture  is  150  m/s  and  the  angle  of  injection  varies  between  0°  and 
60°,  measured  from  the  centerline  of  the  nozzle.  With  this  velocity  and  volume  fraction  a  total 
of  2.01349e-2  Kg/s/rad  of  fuel  are  injected  into  the  domain. 

Full  multiphase  results  are  generated  using  5  droplet  phases  of  sizes  of  60  (tm,  80  pm,  100  pm, 
120  pm,  and  140  pm  with  their  inlet  volume  fractions  being  0.0625,  0.1125,  0.15,  0.1125, 
0.0625,  respectively.  For  the  MUSIG  model,  the  droplet  phase  is  divided  into  10  groups 
using  the  equal  diameter  discretization  with  the  diameter  of  the  smallest  droplet  set  to  55  pm 
and  the  increment  to  10  pm  with  population  fractions  of  0.05,  0.075,  0.1,  0.125,  0.15,  0.15, 
0.125,  0.1,  0.075,  and  0.05,  respectively.  For  the  H-MUSIG  model,  two  variations  are 
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considered.  In  the  first,  two  droplet  phases  are  considered,  each  divided  into  five  group  sizes 
discretized  using  the  equal  diameter  discretization.  In  the  second,  five  droplet  phases  are 
considered,  each  divided  into  two  group  sizes  using  the  equal  diameter  discretization.  The 
diameters  and  population  fractions  of  the  various  groups  are  similar  to  those  used  with 
MUSIG.  The  percentages  of  the  injected  fuel  that  has  evaporated  into  the  gas  field  as 
predicted  by  the  various  methods  are  found  to  be  9.38%,  8.66%,  10.2%,  and  10.1%  for  the 
full  multiphase  method,  the  H-MUSIG  model  with  five  droplet  phases,  the  H-MUSIG  model 
with  two  droplet  phases,  and  the  MUSIG  model  respectively. 

Comparison  of  results  obtained  using  the  various  techniques  are  presented  in  Figures  23 
through  28.  Figure  23  displays  the  velocity  fields  for  some  of  the  droplet  phases.  Figure 
22(a)-22(c)  depicts  the  velocity  vectors  for  droplet  phases  1  (60  pm  in  diameter),  3  (100  pm 
in  diameter),  and  5  (140  pm  in  diameter)  using  the  full  multiphase  approach.  Figures  22(d)- 
(h)  are  for  droplet  velocity  vectors  predicted  by  H-MUSIG  with  Figures  22(d)-22(f) 
displaying  veclocity  vectors  for  the  first,  third,  and  fifth  phases  for  the  droplet  phases  case  and 
Figures  22(g)  and  22(h)  being  for  the  two  droplet  phases  case.  Moreover,  Figure  22(i)  shows 
the  droplet  vector  field  predicted  using  MUSIG.  As  for  the  previous  cases,  larger  droplets 
show  larger  spreading  due  to  their  higher  inertia. 

Comparisons  in  the  form  of  contour  maps  of  the  gas  phase  volume  fraction,  vapor  mass 
fraction,  pressure  field,  and  gas  temperature  fields  generated  by  the  various  methods  are 
depicted  in  Figures  24,  25,  26,  and  27  respectively.  The  general  trend  in  variation  resembles 
the  previous  cases,  i.e.  the  volume  fraction  of  the  particles  decreases  in  the  stream-wise 
direction  (Figure  24),  and  the  mass  fraction  of  the  fuel  vapor  in  the  gas  phase  (Figure  25) 
increases  in  the  streamwise  direction  as  more  kerosene  evaporates.  Pressure  and  temperature 
field  maps  presented  in  Figures  26  and  27  show  similar  profiles.  In  Figure  28,  the  u-velocity, 
gas  temperature,  pressure,  and  vapor  mass  fraction  profiles  across  the  domain  at  x=0.5 
generated  using  the  various  algorithms  are  compared.  As  shown,  the  gas  u-velocity  (Figure 
28(a))  component,  gas  temperature  (Figure  28(b)),  and  pressure  (Figure  28(c))  profiles  are 
nearly  coincident.  Significant  differences  are  notices  in  vapor  mass  fraction  profiles  (Figure 
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28(d))  and  the  increase  in  the  number  of  phases  in  H-MUSIG  has  improved  the  profile  away 
from  the  center  where  it  is  nearly  coincident  with  the  full  multi-phase  profile.  However  in  the 
core  region  the  5  phases  H-MUSIG  solution  predicts  lower  values  than  MUSIG  and  the  two- 
phase  H-MUSIG,  which  are  closer  to  the  full  multi-phase  values  there. 

Closing  Remarks 

Three  numerical  methods  following  a  full  multiphase  approach,  (ii)  a  MUlti-SIze  Group 
(MUSIG)  approach,  and  (iii)  a  Heterogeneous  MUSIG  (H-MUSIG)  approach  for  the 
prediction  of  mixing  and  evaporation  of  liquid  fuel  injected  into  a  stream  of  air  flowing  at  any 
speed  were  developed.  The  numerical  procedures  were  formulated,  following  a  Eulerian 
approach,  within  a  pressure-based  fully  conservative  Finite  Volume  method.  The  k-e  two- 
equation  model  was  used  to  account  for  the  droplet  and  gas  turbulence  with  modifications  to 
account  for  compressibility  at  high  speeds.  The  relative  performance  of  the  three  approaches 
was  assessed  by  solving  for  mixing  and  evaporation  in  three  configurations  involving  droplets 
sprayed  in  the  stream-wise  and  cross-stream  directions  in  subsonic  and  supersonic  streams 
flowing  in  rectangular  and  cylindrical  domains. 

Results,  displayed  in  the  form  of  droplet  velocity  vectors,  contour  plots,  and  axial  profiles 
indicate  that  solutions  obtained  by  the  various  techniques  exhibit  similar  behavior. 
Differences  in  values  are  relatively  small  with  the  largest  being  associated  with  droplet 
volume  fractions  and  vapor  mass  fraction  in  the  gas  phase.  This  is  attributed  to  the  fact  that 
with  MUSIG  and  H-MUSIG  no  droplet  diameter  equation  is  solved  and  the  diameter  of  the 
various  droplet  phases  are  held  constant,  as  opposed  to  the  full  multiphase  approach.  Results 
generated  using  MUSIG  and  H-MUSIG  could  be  improved  through  better  representation  of 
evaporation  in  the  population  balance  equations.  This  will  form  the  subject  of  future 
developments. 
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Figure  Captions 

Figure  1  Schematic  of  the  full  multiphase  approach. 

Figure  2  Schematic  of  the  MUSIG  approach. 

Figure  3  Schematic  of  the  H-MUSIG  approach. 

Figure  4  Physical  domain  for  (a)  streamwise  injection  in  a  rectangular  duct,  (b)  cross- stream 
injection  in  a  rectangular  duct,  and  (c)  streamwise  injection  in  a  cylindrical  duct;  (d) 
an  illustrative  grid. 

Figure  5  Velocity  fields  predicted  by  the  full  multiphase  (a,b,c,  in  increasing  droplet  size), 
the  HMUSIG  (d,  e,  in  increasing  Sauter  diameter)  and  the  MUSIG  (f)  methods  for 
streamwise  injection  in  a  subsonic  flow  field  (Min=0.2)  in  a  rectangular  domain. 

Figure  6  Air  volume  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b), 
and  the  MUSIG  (c)  methods  for  streamwise  injection  in  a  subsonic  flow  field 
(Min=0.2)  in  a  rectangular  domain. 

Figure  7  Vapor  Mass  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b), 
and  the  MUSIG  (c)  methods  for  streamwise  injection  in  a  subsonic  flow  field 
(Min=0.2)  in  a  rectangular  domain. 

Figure  8  Pressure  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b),  and  the 
MUSIG  (c)  methods  for  streamwise  injection  in  a  subsonic  flow  field  (Min=0.2)  in  a 
rectangular  domain. 

Figure  9  Gas  temperature  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b),  and  the 
MUSIG  (c)  methods  for  streamwise  injection  in  a  subsonic  flow  field  (Min=0.2)  in  a 
rectangular  domain. 

Figure  10  Comparison  of  (a)  u-velocity,  (b)  temperature,  (c)  pressure,  and  (d)  vapor  mass 
fraction  profiles  across  the  domain  at  x=0.5m  generated  using  the  full  multi-phase, 
MUSIG,  and  H-MUSIG  methods  (Min=0.2,  Streamwise  injection  in  a  rectangular 
domain). 
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Figure  11  Velocity  fields  predicted  by  the  full  multiphase  (a,b,c,  in  increasing  droplet  size), 
the  HMUSIG  (d,  e,  in  increasing  Sauter  diameter)  and  the  MUSIG  (f)  methods  for 
streamwise  injection  in  a  Supersonic  flow  field  (Min=2)  in  a  rectangular  domain. 

Figure  12  Air  volume  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b), 
and  the  MUSIG  (c)  methods  for  streamwise  injection  in  a  supersonic  flow  field 
(Min=2)  in  a  rectangular  domain. 

Figure  13  Vapor  mass  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b), 
and  the  MUSIG  (c)  methods  for  streamwise  injection  in  a  supersonic  flow  field 
(Min=2)  in  a  rectangular  domain. 

Figure  14  Pressure  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b),  and  the 
MUSIG  (c)  methods  for  streamwise  injection  in  a  supersonic  flow  field  (Min=2)  in  a 
rectangular  domain. 

Figure  15  Gas  temperature  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b),  and 
the  MUSIG  (c)  methods  for  streamwise  injection  in  a  supersonic  flow  field  (Min=2) 
in  a  rectangular  domain. 

Figure  16  Comparison  of  (a)  u- velocity,  (b)  temperature,  (c)  pressure,  and  (d)  vapor  mass 
fraction  profiles  across  the  domain  at  x=0.5m  generated  using  the  full  multi-phase, 
MUSIG,  and  H-MUSIG  methods  (Min=2,  Streamwise  injection  in  a  rectangular 
domain). 

Figure  17  Velocity  fields  predicted  by  the  full  multiphase  (a,b,c,  in  increasing  droplet  size), 
the  HMUSIG  (d,  e,  in  increasing  Sauter  diameter)  and  the  MUSIG  (f)  methods  for 
Cross-stream  injection  in  a  supersonic  flow  field  (Min=2)  in  a  rectangular  domain. 

Figure  18  Air  volume  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b), 
and  the  MUSIG  (c)  methods  for  cross-stream  injection  in  a  supersonic  flow  field 
(Min=2)  in  a  rectangular  domain. 

Figure  19  Vapor  mass  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b), 
and  the  MUSIG  (c)  methods  for  cross-stream  injection  in  a  supersonic  flow  field 
(Min=2)  in  a  rectangular  domain. 
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Figure  20  Pressure  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b),  and  the 
MUSIG  (c)  methods  for  cross-stream  injection  in  a  supersonic  flow  field  (Min=2)  in 
a  rectangular  domain. 

Figure  21  Gas  temperature  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  (b),  and 
the  MUSIG  (c)  methods  for  cross-stream  injection  in  a  supersonic  flow  field 
(Min=2)  in  a  rectangular  domain. 

Figure  22  Comparison  of  (a)  u-velocity,  (b)  temperature,  (c)  pressure,  and  (d)  gas  turbulent 
viscosity  profiles  across  the  domain  at  x=0.5m  generated  using  the  full  multi-phase, 
MUSIG,  and  H-MUSIG  methods  (Min=2,  cross-stream  injection  in  a  rectangular 
domain). 

Figure  23  Velocity  fields  predicted  by  the  full  multiphase  (a,b,c,  in  increasing  droplet  size), 
the  HMUSIG  with  5  droplet  phases  (2  groups  per  phase)  (d,  e,f,  in  increasing  Sauter 
diameter),  the  HMUSIG  with  2  droplet  phases  (5  groups  per  phase)  (g,  h,  in 
increasing  Sauter  diameter), and  the  MUSIG  (i)  methods  for  streamwise  injection  in 
a  supersonic  flow  field  (Min=2)  in  a  cylindrical  domain. 

Figure  24  Volume  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  with  five 
droplet  phases  (2  groups  per  phase)(b),  the  HMUSIG  with  tow  droplet  phases  (5 
groups  per  phase)(c),  and  the  MUSIG  (d)  methods  for  streamwise  injection  in  a 
supersonic  flow  field  (Min=2)  in  a  cylindrical  domain. 

Figure  25  Vapor  mass  fraction  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  with 
five  droplet  phases  (2  groups  per  phase)(b),  the  HMUSIG  with  tow  droplet  phases 
(5  groups  per  phase)(c),  and  the  MUSIG  (d)  methods  for  streamwise  injection  in  a 
supersonic  flow  field  (Min=2)  in  a  cylindrical  domain. 

Figure  26  Pressure  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  with  five  droplet 
phases  (2  groups  per  phase)(b),  the  HMUSIG  with  tow  droplet  phases  (5  groups  per 
phase)(c),  and  the  MUSIG  (d)  methods  for  streamwise  injection  in  a  supersonic 
flow  field  (Min=2)  in  a  cylindrical  domain. 
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Figure  27  Gas  temperature  fields  predicted  by  the  full  multiphase  (a),  the  HMUSIG  with  five 
droplet  phases  (2  groups  per  phase)(b),  the  HMUSIG  with  tow  droplet  phases  (5 
groups  per  phase)(c),  and  the  MUSIG  (d)  methods  for  streamwise  injection  in  a 
supersonic  flow  field  (Min=2)  in  a  cylindrical  domain. 

Figure  28  Comparison  of  (a)  u-velocity,  (b)  temperature,  (c)  pressure,  and  (d)  vapor  mass 
fraction  profiles  across  the  domain  at  x=0.5m  generated  using  the  full  multi-phase, 
MUSIG,  and  H-MUSIG  methods  (Min=2,  Streamwise  injection  in  an  Axi- 
Symmetric  domain). 


